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iH ■ Abstract 

We investigate the mechanism in the tunnehng dynamics of open ultracold few-boson systems 
by numerically solving the time-dependent few-boson Schrodinger equation exactly. By starting 
from a weakly repulsive, initially coherent two-boson system we demonstrate that the decay dy- 
namics incorporate fragmentation. The wavefunction of the tunneling state exhibits a pronounced 
dynamically-stable pattern which we explain by an analytical model. By studying more bosons 
and stronger interactions we arrive to the conclusion that the decay by tunneling is not a coherent 
process and exhibits a wealth of phenomena depending on the interaction between the particles. 
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One of the most basic and also paradigmatic phenomena in quantum mechanics is the 
tunnehng process in open systems on which ample studies have been reported (for a recent 
l|). Much attention has also been devoted to the tunneling of particles between 



book, see 



m 



potential wells (see, e.g., |2l, |3|), but in the present work we concentrate on open systems. 



where the tunneling is from a wel 



m 



through a potential barrier into open space. Since Bose- 



Einstein condensates (BECs) j^, |5| have been realized, this unique state of matter has been 



used to model a wealth of physical phenomena in a controllable fashion [6[. BECs were made 
one-dimensional, see, e.g., [Tf, and there is progress towards the realization of ultracold few- 
boson systems 8|. Such systems provide ideal real objects to investigate the mechanism 
behind the tunneling process in interacting systems and to arrive at a deeper understanding 
of the quantum dynamics of tunneling in general. 

When studying tunneling the calculation is often reduced to that of a single particle. In 
the case of BECs the time-dependent Gross-Pitaevskii (TDGP) theory, which is known to 
be successful in many situations, naturally describes tunneling of the interacting bosons as 
the tunneling of a single effective boson. There are very few studies available which discuss 
tunneling in open systems beyond the Gross-Pitaevskii (GP) theory. For infinitely strong 
interacting bosons (Tonks-Girardeau gases), Del Campo et al. |9| analytically treated the 



case of decay by tunneling through a Dirac-(5 barrier. In [10| the decay rates of a BEG in 
a double-barrier potential were computed by employing complex scaling to the GP and the 
best mean field methods. This investigation provided first indication that tunneling of an 
initially coherent system can lead to fragmentation, which is a known phenomenon of BECs 
confined to an external potential Uj. This finding has recently been substantiated in [12] 
where it is demonstrated that the TDGP theory fails to correctly describe the tunneling 
of open few-boson systems even for weakly interacting particles. Interestingly, this failure 
is more severe for four than for two particles [12b]. Clearly, the mechanism governing the 
physics of tunneling in open few-boson systems involves a dynamical loss of coherence of the 
wavefunction, a phenomenon which cannot be described by GP theory. 

In this work we solve the time-dependent Schrodinger equation (TDSE) numerically ex- 
actly, and unravel the mechanism of time-dependent tunneling in one dimension. We demon- 
strate fragmentation even for weak interactions. A particularly exciting stationary pattern 
of the tunneled fraction emerges which is analyzed by a simple analytical model. We start 
with two bosons and show that the dynamics in the case of weak interactions exhibit dy- 



namical fragmentation: the density forms a dynamically-stable interference-like pattern, the 
momentum distribution exhibits two prominent peaks, and the largest occupation number 
decreases, i.e., the system becomes less condensed (coherent). All these features cannot be 
covered by TDGP theory. Increasing the interparticle repulsion, we verify that the density 
mimics fermionic behavior, expected in the limit of very strong interactions, see work by 
Girardeau 13|. Next, considering four particles we see that the dynamically- stable struc- 
ture found for two bosons disbands and becomes time-dependent, but the phenomenon of 
dynamical fragmentation remains. 

We tackle the TDSE with the full many-body Hamiltonian, H = Xlfci^l^j) + 
X] ■<k ^{^j ~ ^k) in dimensionless units. Here the one-body Hamiltonian is given by 
h{x) = T{x) + V{x). T is the kinetic energy operator T(x) = — 19^ and V{x) is the 
potential energy operator, shown by the dashed and solid black lines in Fig. [TJ\ (for an an- 
alytical form see Ref. [12]). To model the commonly used Dirac-(5 two-body interaction we 
take W{xj—Xk) = \oSa{xj—Xk) where 6a{x) is a normalized Gaussian, of width a = 0.05. Aq 
is the interaction strength. To emulate an open, infinite system we use a complex absorbing 
potential (CAP) |l4l. Il5l|. which we add to the one-body Hamiltonian h(x). It absorbs the 
wavefunction from some point Xs on and is adjusted such that it is essentially reflection 
free and perfectly absorbing (see [14]). We use the multi-configurational time- dependent 



Hartree (MCTDH) ly, [iTJ] method in its parallel version |18| to solve the few-boson TDSE 



numerically exactly. This algorithm has already been successfully applied to the dynamics 



of few-boson systems (see, e.g., |l2Lll9|). The parameters for convergence of the calculation 
with respect to the CAP, number of grid points and the number of orbitals are the same as 
in Ref. [12]. In the present work, we enlarge the simulated region to x G [—6,226], and use 
2048 grid points. The CAP starts at xs = 196. 

The most general quantity to describe a quantum mechanical system is the A^-body 
density which is a function of as many variables as the system's constituents, A^, plus one 
for the time. The simplest interacting system to consider is that of two weakly interacting 
bosons. Moreover, in this case we are able to visualize the full density, |\I'(xi,X2; t)p, as a 
3D-plot at a given time t = t'. Thus, we prepare the ground-state of the weakly interacting 
two-boson system with A = Ao(A^— 1) = 0.3 in a parabolic potential (see Fig.[T]A). This initial 
state is almost completely coherent, its largest occupation number pi{t = 0) is 99.85%. Then 
we abruptly switch the parabolic trap to the open potential. Consequently, the wavepacket 



is decaying by tunneling through the barrier. In Figs. [T] B-D we plot the density of the 
exact solution for the two-boson case with interparticle interaction strength A = 0.3 at 
times t = 0; 300; 400, respectively. In Fig. [IP the TDGP density at t = 400 is depicted 
for comparison. First, we notice that a dynamically stable shape is formed and persists 
for long time on the scale of the half-life of the system. We refer to this shape of the 
density as the lobe-structure. Second, the TDGP density, contrasting the exact solution, 
does not display such a lobe-structure. Several questions arise: Is the tunneling process 
a coherent process? Does the initial almost fully coherent, weakly interacting system stay 
coherent throughout its time evolution? Surprisingly the answers are negative. The fact 
that the computed exact density differs drastically from the TDGP density allows us to 
interpret this difference as a loss of coherence. We recall that in the TDGP theory the 
wavefunction stays completely coherent at all times, i.e. pi{t) = 100%, Vt. Thus we have 
to ask: What is the mechanism behind the system's loss of its coherence and what makes 
the system form a dynamically stable shape in the density? Does this shape persist for 
more particles or stronger interactions? To answer the former question we will investigate 
the momentum distribution of the tunneling dynamics and to answer the latter we will 
subsequently look at the density of a strongly interacting two-boson system and at the 
two-body- density p^'^\x'i = xi, x'2 = X2,xi,X2;t) = p^'^\xi,X2; t) 20| of a four-boson system. 



To get a deeper insight on the many-body quantum dynamics we investigate the momen- 
tum distribution of the considered system. In Fig. [2]we compare the momentum distributions 
of the exact solution and the TDGP for the times t = 0, 300, respectively. We see that the 
momentum distributions of the initial guesses are similar and Gaussian-shaped, as it is ex- 
pected for weakly interacting bosons in a parabolic trap. For the time-evolution, in turn, 
we get a persistingly prominent two-peak structure on a Gaussian background for the exact 
solution and a persistingly prominent single-peak structure for the TDGP dynamics. We 
attribute the peak-structure to the part of the density running away from the well. Thus, 
we can conclude that the single momentum dominated picture of the dynamics is caused 
by the constraint that the system stays in a single quantum mechanical state, i.e., is fully 
coherent. But the exact density shows two peaks in the momentum distribution. Hence, 
we can suppose that the system occupies two states, exhibiting two distinct momenta. A 
bosonic system occupying two quantum mechanical states is fragmented - this allows us 
to state that the loss of the initial coherence is caused by the tendency of the system to 



fragment in the course of the tunnehng dynamics outside the barrier. 

We observe the lobe-structure in the density for a whole range of interactions, < Aq ^ 2. 
We have also found that the distance of the peaks increases with the increase of repulsion, 
whereas the length scale of the lobe-structure decreases. These observations motivated us to 
formulate an analytical model for the wavefunction. We assume that in the outside region, 
i.e., at X > xc, the wavefunction is in general constituted by plane waves. The repulsion 
forces the bosons to occupy two plane waves, i.e., (f)j{x,t) = (1 — Qxc)'^in,jit)(f)inj{x) + 
Qxc^-jifye^^^^^'ii = 1, 2. The 0i„j- are the ground and first excited state of a single particle in 
a harmonic potential. The coefficients aj„j(t), aj{t) describe the time-dependent occupations 
of the inner parts, (j)in.j{x,t), and the plane waves, respectively. The Heavyside function, 
Q^^ = 6(x — xc), starts at the barrier, xc = 2.25, and the kj are the positions of the two 
peaks occurring in the momentum distribution of the tunneling dynamics. Assuming further 
that the coefficients Oj (t) are real we can construct the full density from these orbitals: 



p^^>{xi,X2;t) = {l-eccc)pin{xuX2;t) 



+ Qxc { Mt) COS [{k2 - ki) xi] COS [{k2 - ki) X2] (1) 



+ B{t) cos [{k2 - ki) (xi + X2)] + C{t) \. 

Here the coefficients A{t), B(t) and C(t) collect the dependency on the coefficients aj(t), 
and p\J is a symmetrized product of the (pinj- In Fig. [3] we plot p^'^\xi,X2;t) and take 
{k2 — ki) = 0.15 from Fig. |2] for the distance of the momentum peaks and A(t) = B{t), 
C{t) = 0. We note here that this choice of A, B and C corresponds to a fully fragmented 
wavefunction outside the well. Comparing Fig. [3JA. with Fig. [33 and C, we find that our 
model reproduces very well all key features found. The lobe-structure is stable in time for a 
range of interactions. Another nice feature of the model is that it also reproduces the TDGP 
results as a limiting case. Indeed, if we consider a single 0j(x, t) [take ki = ^2 in Eq. ([1])], we 
obtain the single value of A{t) + B(t) + C{t) in the outside region. So, the non-oscillatory, 
spatially-independent behavior of the TDGP density in Fig. [T]E is reproduced. We conclude 
that the mechanism behind the loss of coherence of the system is that the tunneled density 
dynamically fragments and occupies mainly the permanents formed by two plane waves 
given in Eq. ([1]). Thus, the lobe-structure is a self-interference pattern. Clearly this process 
cannot be described at all by the TDGP theory. 



The occurring fragmentation can also be observed by looking at the time evolution of the 
i = 1, ...,M renormalized natural occupation numbers Pi{t). In Fig. [3l3 we plot their time 
evolution. The fragmentation appears as an increase of P2{t) to about 10% accompanied by 
a decrease of pi(i) to around 90%. The other occupation numbers, Pi>2{t), stay very small 
at all times. This describes the loss of coherence in the tunneling dynamics from another 
viewpoint. Let us analyze the dynamics of the occupation numbers a little closer: The 
initial time evolution is accompanied by oscillations of {pi{t);i = 1,2}. From our model 
consideration we can conclude that the relative occupancy, ai(t) / a2(t) , of the plane waves 
constituting the outside fraction of the wavefunction does not change much in time - that is 
why we find a dynamically-stable density. Therefore, the mentioned oscillations are caused 
by the density localized in the well. We term the oscillations beating. We find that the 
frequency of the beating, u, is related to momentum by u; = (fci — ^2)^/2, where |fci — ^2] 
is the spacing of the two peaks in the momentum distribution in Fig. [2l Thus, the beating 
and the appearance of the two peaks in the momentum distribution are intimately related. 

We analyze the dynamics of two more strongly repelling bosons, to see whether and how 
the above discussed mechanism of tunneling dynamics is altered. When we increase the inter- 
boson repulsion further (A > 2), the dynamically-stable structure dissolves and becomes 
time-dependent. According to the reasoning of our model this happens when more than 
two states outside the well become populated. Obviously, the process of loss of coherence 
intensifies in this case as the wavefunction occupying more than two states is less coherent 
than the one occupying two. When we increase the interaction strength to the limit of very 
strong repulsion, e.g., A > 60, the density starts to mimic the one in the fermionization limit 



see 



13(1). In Fig. |3p we plot the computed density for A = 60 and A^ = 2 at t = 150. It is 
clearly seen that the density's diagonal becomes 0, i.e., \l/(a;i = X2',t) = 0. The off-diagonal 
shows strong dependence on time. 

Let us see now how the found mechanism behind the tunneling process develops when 
increasing the particle number. We investigate a four-boson system with an equal mean 
field parameter A = Ao(iV — 1) = 0.3. We recall that all systems with the same A reveal 
identical TDGP dynamics, irrespective to changes in the number of bosons, N. Now, i.e., for 
A^ = 4, we consider the diagonal part of the reduced two-body density matrix, p*^^)(xi, X2; t), 
in order to be able to compare with N = 2, where the density \'${xi,X2]t)\'^ is equal to 
p(^)(a;i,X2; t). In Fig. [HP we plot p^'^^ ioi N = A at t = 200. The shape of the density 



exhibits time-dependency in the outside region. This indicates that adding more bosons 
to the system also forces it to occupy more than two states during the dynamics. We find 
that the dynamically-stable self-interference pattern is not persisting as the process of loss 
of coherence intensifies when we increase A^. This explains why the survival probability 
computed on the many-body level deviates from the TDGP results even more for four than 
for two bosons [12b]. 

Let us conclude. We have identified the mechanism governing the tunneling dynamics in 
open ultracold bosonic systems based on an exact numerical solution of the TDSE. The re- 
sults were compared to the TDGP approximation, which is widely and successfully used, but 
is not applicable to the scenario studied in this work. We explain the considerable deviations 
of the exact solution from the TDGP dynamics by the fragmentation of the wavefunction, 
occupying more states in the outside region. For the case of two weakly interacting bosons 
we observe a lobe-structure which is dynamically-stable. A simple analytical model of a 
two-fold fragmented state formed by two plane waves with different momenta explains this 
lobe-structure. 

For stronger interactions for two bosons and for a larger number of bosons the 
dynamically-stable structure breaks and becomes time-dependent, indicating thereby that 
the loss of coherence intensifies. The tunneling of a few interacting bosons is not a coherent 
process and the interference effects taking place in the outgoing fraction of the wavepacket 
are in this case a signature of fragmentation. 
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FIG. 1: (Color online) Tunneling dynamics of weakly interacting bosons: A = 0.3 and N = 2 
bosons. (A): The initial parabolic trap is depicted by the black-dashed line. The potential used 
for the propagations is plotted as solid black line. The density at times t = 0, 300, 400 is indicated 
by the green, blue and magenta line, respectively. The energy per particle is plotted as a turquoise 
horizontal line. The inset shows an enlargement of the density in the outside region. (B),(C),(D): 
The wavefunction (square) evolves from its initial Gaussian shape to a lobe- structure which is 
persistent throughout the dynamics. (See text for details.) (E): For comparison, the TDGP shows 
minor oscillations and no lobe- structure. All quantities shown are dimensionless. 
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FIG. 2: (Color online) The momentum distribution for the times t = for the exact solution 
(black line) and t = 300 for the exact solution (red line) and the TDGP solution (green line). The 
initial guesses are coherent, i.e., at t = the exact and the TDGP momentum distribution cannot 
be distinguished on the scale of the plot. The TDGP dynamics have a single-peaked momentum 
distribution. The exact momentum distribution reveals a two-peak structure. The intensity of the 
peak in the TDGP momentum distribution is roughly the sum of the intensities of the two exact 
peaks. All quantities shown are dimensionless. 
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FIG. 3: (Color online) Self-interference and fragmentation is the mechanism governing the tunnel- 
ing dynamics. (A): The model introduced in Eq. ([T]) for the wavefunction (square) is plotted for a 
fully fragmented outer wavefunction. Note the similarity to the dynamically-stable structure de- 
picted in Figs.dl^C) and (D). (B): The time-evolution of the occupation numbers {pi{t);i = 1,2,3} 
(normalized to 1 Vt) for N = 2 weakly interacting bosons indicate the dynamical fragmentation of 
the system. (C): The wavefunction (square) |^(xi,X2;i = 150)|^ for a rather strong interaction of 
the 2 bosons. The diagonal, xi = X2, tends to 0. We show only one specific time, t = 150, here, 
but in general the off-diagonal part shows short-range oscillations and strong time-dependency. 
(D): The two-body density, p^'^'{xi,X2',t = 200), of A^ = 4 weakly interacting bosons. Although 
we depict only t = 200, let us state that p^"^' is time-dependent and forms no dynamically-stable 
pattern. See text for discussion. All quantities shown are dimensionless. 
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